# geometry_optimization.py
# Geometry optimization of camphorquinone
# Method: B3LYP/6-31G(d)
# Solvent: PCM (epsilon = 4.0)

from pyscf import gto, dft
from pyscf.solvent import pcm
from pyscf.geomopt.geometric_solver import optimize

# Initial geometry for camphorquinone (C10H14O2)
# Generated from SMILES: CC1(C)C2CCC1C(=O)C2=O
# Pre-optimized with UFF force field
mol = gto.Mole()
mol.atom = '''
C    -0.825830    1.197110    0.284570
C     0.583850    0.704310   -0.113320
C     1.378790    1.952940   -0.530900
C     0.460290    3.168390   -0.355700
C    -0.896510    2.717100   -0.009530
C    -1.730930    2.436250    1.249270
C    -1.659150    0.207610   -0.214040
O    -2.799260   -0.090010    0.138970
C     0.660540    0.015690   -1.483840
O     0.831490   -1.183600   -1.650210
H     0.601980   -0.010090    0.715980
H     2.276390    2.120670    0.085880
H     1.702420    1.835760   -1.578830
H     0.829910    3.926130    0.353660
H     0.426430    3.678530   -1.330660
H    -1.510490    3.484710   -0.517210
H    -1.220110    1.711880    1.909730
H    -1.950050    3.364540    1.801540
H    -2.692540    1.981180    0.998580
H    -1.082900   -0.311590   -1.020110
H     0.551720    0.636210   -2.367630
'''
mol.basis = '6-31g(d)'
mol.charge = 0
mol.spin = 0
mol.build()

print("="*70)
print("CAMPHORQUINONE GEOMETRY OPTIMIZATION")
print("Method: B3LYP/6-31G(d)")
print("Solvent: PCM (ε = 4.0)")
print("="*70)
print()

# DFT calculation with B3LYP
mf = dft.RKS(mol)
mf.xc = 'B3LYP'
mf.verbose = 4

# Add PCM solvent
mf = pcm.PCM(mf)
mf.with_solvent.eps = 4.0
mf.with_solvent.verbose = 4

print("Starting geometry optimization...")
print()

# Geometry optimization with tight convergence
mol_eq = optimize(mf, maxsteps=100)

print()
print("="*70)
print("OPTIMIZATION CONVERGED")
print("="*70)
print()

# Get final SCF energy
mf_final = dft.RKS(mol_eq)
mf_final.xc = 'B3LYP'
mf_final = pcm.PCM(mf_final)
mf_final.with_solvent.eps = 4.0
E_final = mf_final.kernel()

print(f"Final SCF Energy: {E_final:.10f} Hartree")
print()

# CRITICAL: Write optimized geometry to XYZ file
with open("camphorquinone_optimized.xyz", "w") as f:
    f.write(f"{mol_eq.natm}\n")
    f.write("Camphorquinone optimized geometry | B3LYP/6-31G(d)/PCM(ε=4.0)\n")
    for atom in mol_eq.atom:
        symbol, coord = atom[0], atom[1]
        f.write(f"{symbol} {coord[0]: .8f} {coord[1]: .8f} {coord[2]: .8f}\n")

print("="*70)
print("OUTPUT FILES GENERATED")
print("="*70)
print("✓ camphorquinone_optimized.xyz")
print()
print("Next step: Run frequency_analysis.py to verify this is a true minimum")
print("="*70)
